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Resonant tunneling and intrinsic bistability in twisted graphene structures 

J. F. Rodriguez-Nieva^, M. S. Dresselhaus^’^, L. S. Levitov^ 

^Department of Physies, Massaehusetts Institute of Teehnology, Cambridge, MA 02139, USA and 
^Department of Eleetrieal Engineering and Computer Seienee, 

Massaehusetts Institute of Teehnology, Cambridge, MA 02139, USA 

We predict that vertical transport in heterostructures formed by twisted graphene layers can 
exhibit a unique bistability mechanism. Intrinsically bistable I-V characteristics arise from resonant 
tunneling and interlayer charge coupling, enabling multiple stable states in the sequential tunneling 
regime. We consider a simple trilayer architecture, with the outer layers acting as the source 
and drain and the middle layer floating. Under bias, the middle layer can be either resonant or 
non-resonant with the source and drain layers. The bistability is controlled by geometric device 
parameters easily tunable in experiments. The nanoscale architecture can enable uniquely fast 
switching times. 


I. INTRODUCTION 

Nanoscale systems that can switch between distinct 
macroscopic states upon variation of some control param¬ 
eter are in high demand in diverse areas of nanoscience 
research. Bistable electronic systems which exhibit fast 
switching are of interest for applications such as low- 
power memory and logic [1] . Recently, new realizations of 
intrinsically bistable system have been discovered, both 
in graphene [m and in other systems UMl- In particu¬ 
lar, van der Waals heterostructures comprising graphene 
layers sandwiched between insulating hexagonal boron- 
nitride (hBN) layers afford electronic environments with 
tailored band structures and transport characteristics 
[To]. It was demonstrated that introducing a twist be¬ 
tween adjacent graphene layers in such heterostructures 
can result in a resonant behavior of the tunneling cur¬ 
rent and non-monotonic I-V characteristics m- It is 
therefore tempting to exploit twisted graphene multi¬ 
layer structures as a platform for bistable and hysteretic 
nanoscale systems. 

Here we predict intrinsic bistability and hysteretic I-V 
characteristics for vertical transport in heterostructures 
formed by graphene monolayers separated by hBN bar¬ 
riers, in a twisted arrangement similar to that described 
in Ref. m- Essential for our bistability mechanism are 
resonances originating from momentum-conserving tun¬ 
neling between linearly dispersing Dirac bands m and 
occurring when the bands are aligned m (see Figl^,c). 
Bistability arises due to current-induced charge accumu¬ 
lation producing an interlayer bias that tunes the inter¬ 
band tunneling in and out of resonance. 

Below we focus on the simplest case of a two-step se¬ 
quential tunneling in a device comprising three graphene 
monolayers. Such trilayer architecture, pictured in 
FigHt, with the top and bottom layers acting as a source 
and drain and the middle layer electrically decoupled 
(floating), is similar to previously studied double-barrier 
quantum-well (QW) structures [13]. However, our bista¬ 
bility mechanism, originating from resonant tunneling 
between Dirac bands in graphene layers, is distinct from 


that in the QW structures [13]. In our case, multiple 
stable states arise because the decoupled layer can, for 
a fixed external bias, be either in a resonant (low resis¬ 
tance) or a non-resonant (high resistance) state. This 
behavior is illustrated in Figj^. 

The bistability is governed by geometric parameters 
- the twist angle 0 and the interlayer distances dij - 
which are easily tunable in experiments. The twist angle 



FIG. 1. (a) Trilayer graphene heterostructure schematics, 

with layers labeled 1 to 3. Here Ej and dij are the interlayer 
currents and distances, (b) Band structure of the twisted 
graphene layers l(blue) and 2(red). The twist angle 0 de¬ 
fines a characteristic energy A [Eq.Q] and three superlattice 
wave vectors qA,n,c [Eq. (c) Bistable I-V characteris¬ 

tics. The resonant and non-resonant bistable states are illus¬ 
trated in the top left inset (details are discussed in Figj^. 
The procedure for hnding bistable solutions is illustrated in 
the bottom right inset [see Eq. and accompanying discus¬ 
sion] . 
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controls the Dirac cones’ displacement in the two layers 
and the energy at which the cones intersect (see 

Iq^l = (87r/3ao)sin(0/2), A = (1) 

where v-p 10® m/s is the carrier velocity and oq ~ 
2.46 A is the graphene lattice constant. The distances 
dij^ marked in Figj^, determine the interlayer tunnel 
conductance values Gij ^ where A is the WKB 

length governing the tunneling amplitude dependence on 
barrier width. In what follows we will use the conduc¬ 
tance ratio 

Z = G 12 /G 23 - e2(‘^23-di2)/A (2) 

where Gij denotes the conductance between the corre¬ 
sponding layers. 

The quantities 0 and dij can be controlled with a large 
degree of precision. The twist angle 0 can be tuned within 
^ 1° during fabrication m , whereas dij can be varied by 
adding monolayers of dielectric materials, such as hBN 
or M 0 S 2 . Since typical values A = hl{2meW)^^‘^ ^ 2 A, 
estimated for the tunneling barrier height IT ^ 1 eV and 
the effective electron mass rrie ^ 10“^^ kg, are compara¬ 
ble to the hBN or M 0 S 2 monolayer thickness, variation 
in dij results in a fairly gradual change of Z. 

One appealing aspect of this system is the short inter¬ 
layer transport length of a nanometer scale, which can 
allow high operation speeds and fast switching times. 
This is evident from an estimate for the RG time, 
'Trc = njAiigd ^ 100 ns, where r ^ 1 is the dielec¬ 
tric constant, d ^ 1 nm is the interlayer separation, and 
g ^ is the interlayer conductance per 

unit area. The combination of geometric tunability and 
small transport lengths is not present in previously stud¬ 
ied graphene-based bistable systems, such as graphene 
flash memories 0131 or graphene resistive memories HE]. 
Small thicknesses can also enable large packing densities. 

The steep electronic dispersion in graphene makes the 
bistable state properties distinct from those in QW sys¬ 
tems. In our case, the bistability is controlled by the res¬ 
onances arising due to band alignment. The correspond¬ 
ing bias value, which scales as a power law of the energy 
A given in Eq. 0, can be as large as SV ^ 100-500 mV 
(see discussion below). In QW systems, instead, the bias 
range where bistability occurs is mainly controlled by the 
amount of charge ngw that can be stored in a quantum 
well, 6V « euQw/G^ where G is the interlayer capac¬ 
itance [T4|. Typical carrier densities in the “charged” 
and “uncharged” states of a bistable QW system, as¬ 
sessed by magnetic oscillation measurements [T5|, are on 
the order of ngw ^ 10^^/cm^ and ngw ^ 0, respectively. 
These carrier densities yield typical values 6V ^ 50 mV 
in double-barrier quantum wells with a width of tens of 
nanometers {G ^ 0.1-1 mF). Such values can be as much 
as an order of magnitude smaller than the above estimate 
predicts for the graphene case. 


II. SEQUENTIAL TUNNELING MODEL 

Vertical transport in our trilayer architecture can be 
described by a simple sequential model. The model va¬ 
lidity relies on the interlayer tunnel coupling being weak 
such that the inter-layer charge transfer is slow compared 
to the intra-layer electron relaxation. Indeed, the val¬ 
ues trc^ estimated above, are much longer than typical 
thermalization times in graphene, rth ^ 10 ps [T6|. The 
RG times, however, are sufficiently fast to be competitive 
with the speeds of existing switching devices [1] . 

The interlayer transport mechanism is mainly gov¬ 
erned by the twist angle 6>, which defines the K-point 
displacement between graphene lattices in adjacent 
layers, and the interlayer bias. Under bias, the value 
|qy^| given in Eq. 0 determines the range of momenta 
and energies for which momentum-conserving tunneling 
is allowed. Large values of |q^| hinders resonant tunnel¬ 
ing given that phonon and defect scattering are necessary 
to supply the large momentum mismatch between layers. 
Momentum-nonconserving transport can also occur if the 
top/bottom layer is made of a different material so that 
there is a large mismatch between unit cells with respect 
to that of graphene. For small |q^|, on the other hand, 
momentum conserving tunneling is possible for moder¬ 
ately small values of bias. 

In our two-step sequential tunneling model, we treat 
transport between layers 1 and 2 as momentum- 
conserving. The second step, between layers 2 and 3, is 
assumed to be momentum-nonconserving and described 
by Ohm’s law. The latter assumption allows us to sim¬ 
plify our discussion and focus on the essential aspects of 
bistability. In addition, we also assume that the contact 
resistances are sufficiently small so that all the potential 
drop occurs predominantly between the graphene layers. 

Turning to a systematic development of the model, the 
low energy Hamiltonian R describing coherent transport 
between a pair of twisted graphene monolayers has con¬ 
tributions R = Ri -f R 2 + 7i 2- Here Hi ,2 are the free- 
particle terms describing massless Dirac particles in each 
graphene layer, and T 12 describes the interlayer tunnel 
coupling HMS]. The free particle terms are 

^1 = E V’fkl^^Fcr ■ (k + q2i/2) - 

^ (3) 

”^2 = E “ qA/2) - M2]^/’2.k, 

k 

where /ii ,2 are the Fermi energies measured relative to 
the Dirac point. For a small twist angle 6>, the large- 
wavenumber processes that couple different valleys can 
be neglected. In this case, it is sufficient to account for 
a single Dirac cone in each layer, see Eq. 0. We adopt 
this approximation below. 

The tunneling coupling can be modeled as a local, pe- 
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FIG. 2. Twisted graphene layers form an hexagonal superlat- 
tice with reciprocal superlattice vectors Gi and G 2 [l7]. The 
momentum conserving tunneling coupling has the periodicity 
of the superlattice and can be decomposed into Fourier com¬ 
ponents G = nGi + mG 2 , with n^m being integers. For a 
small twist angle 0, tunneling is dominated by the smallest 
wavevectors qa, qs = Qa —Gi and qc = qA —G 2 , see Eq. 

riodic function of position HI]: 

Ti2 = ^'0pk'^G'^2,k+G + H.C.. (4) 

k,G 

The periodicity of the interlayer coupling, quantified by 
the G wavevectors, is determined by the hexagonal su¬ 
perlattice unit cell that is formed by the twisted graphene 
layers, see Figj^ For small only the longest wavelength 
contributions are relevant for tunneling. Referred from 
the Dirac point of layer 1, such long wavelength compo¬ 
nents are given by qA, Qb = Qa — and qc = Qa — G 2 
(see Fig§, where Gi ^2 are the reciprocal vectors of the 
superlattice Brillouin zone, which is smaller than the 
graphene Brillouin zone by a factor ^ sin^(^). While the 
higher-q harmonics of the interlayer hopping potential 
spatial modulation also contribute to tunneling, it can 
be shown that their contributions vanish rapidly on the 
reciprocal lattice vector Gi ,2 scale [IZlIIH]. This leads to 
the tunneling Hamiltonian 


e.g. in highly-oriented hBN-graphene structures, which 
have a small lattice mismatch of about 1.8% (a detailed 
discussion of these effects can be found in Ref. [20] ). Such 
effects, if present, would alter the values qA(B,c) but oth¬ 
erwise not change our discussion in an essential way. 

Under an interlayer bias potential V 12 , the tunneling 
current /12 is 


/l2 = 


eN 


X A2,5'(k + qj, w) [/i(‘^) - / 2 (w)], 


(6) 


where s (s') refers to the electron (-b) and hole (—) bands 
of layer 1 (2), and = 4 is the spin and valley degener¬ 
acy. The functions fi{uj) = ^ 1 ] the Fermi 

distribution functions for each layer, with (3 = 1 /k^T 
being the inverse thermal energy and /i^ being the Fermi 
energies. The function Ai^g is the spectral function of 
layer i and band s. The energy for the quantities in layer 
2 is offset by d) = cj + e4>i2 due to the built-up interlayer 
electrostatic potential 4>i2 [see Eq. (§] between layers 1 
and 2. Because of capacitance effects, the interlayer elec¬ 
trostatic and chemical potentials are related by 


eVi2 = Ml - M2- e$i2, (7) 


where and 4>i2 are implicit functions of V 12 . The quan¬ 
tity Tj^ in Eq.Q denotes 
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(k) = (k,s, l|Tj|k + qj',s',2), 


( 8 ) 


where |k, s, i) are the two-component eigenvectors of Hi ,2 
in Eq.(|^ and 6>k is the k- vector polar angle. 

The bistability can now be described by combining re¬ 
lations @ and 0 as follows. In a steady state, there 
is zero net flow of carriers into the middle layer. There¬ 
fore, when the external bias V = V 12 + V 23 between top 
and bottom layers is fixed, the equilibrium current I is 
obtained by solving for V 12 from the non-linear equation 


/(U)=/i2(Ui2)=/23(^-V^12). (9) 


7i 2 — ^l.k "^3 V^2,k+q^- + H.C. (5) 

j=A,B,C k 

comprised of only three Fourier components. In this ex¬ 
pression for 7 i 2, the k vectors are relative to the Dirac 
point of each layer, i.e. k — qA/2 —k in layer 1 and 
k + qA/2 ^ k in layer 2. 

Parenthetically, the lattice of the dielectric material 
separating the graphene layers can produce slowly vary¬ 
ing spatial modulation of the tunneling transition am¬ 
plitude T in Eq.([^, giving rise to the effects resembling 
those due to a twist angle 0. This would be the case when 
the dielectric and graphene are nearly lattice-matched as 


This procedure to obtain the I-V response is shown 
graphically in the inset of Fig{^. The straight line de¬ 
scribes transport between layers 2 and 3 which is assumed 
to follow Ohm’s law, I 23 = G 23 V 23 , where G 23 and V 23 
are the interlayer conductance and interlayer bias poten¬ 
tial between layers 2 and 3, respectively. 

III. ELECTROSTATIC FEEDBACK 

In order to include the electrostatic feedback, Eq. 0 
needs to be complemented with further electrostatic con¬ 
siderations that relate the variables Uj, and /i^. It 
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is important to note that all variables can be determined 
once the carrier densities in each layer, ni, n 2 and ns, 
are known. Indeed, assuming that there is no external 
gate, the neutrality condition relates the charge densities 
in the different regions of the device as 

+ ^2 + ^3 = 0 . ( 10 ) 

Furthermore, the application of an external bias potential 
V fixes the Fermi level difference between layer 1 and 
layer 3 as 


47re^ 

eV = Jill - jLi 3 -i -(nidis + n2d23). (H) 


Here dij is the interlayer distance between layer i and 
layer j, k is the dielectric constant of the barrier mate¬ 


rial, and iJ^i = sgn(n^) hv ^. In Eq.(ll), we implicitly 
assume that all layers are undoped at H = 0. Equa¬ 
tions @-0 then form a closed set of equations from 
which ni, n 2 and ns can be obtained. The remaining 
variables, Vij and are functions of n^. In particu¬ 
lar, the electrostatic potentials are 4>i2 = —47re^di2ni//^ 
and 4>2 s = 47re^(ni(iis + '^ 2 <^ 23)//^5 whereas the inter¬ 
layer bias potentials are V 12 = — (/i 2 + ^ 12 ), and 

^23 = (/^2 + ^12) — (M3 + ^23)- 

Eor simplicity, here we fix the Eermi energies in Eq. ^ 
to a constant value fii = /i. This is equivalent to turn¬ 
ing off capacitance effects. In this case, Vij = 

(see EigJ^. This approximation is valid in the regime 
Ae^dij/S^l ~ 15 • dij[nm]A[eV]//^ > 1. In this 
regime, minimal changes in carrier concentration induce 
large interlayer electrostatic potentials. The more re¬ 
alistic scenario which includes quantum capacitance ef¬ 
fects ED, such that /ii ^2 vary with V 12 , is here consid¬ 
ered in Appendix A. However, this more realistic picture 
only introduces small corrections to the tunneling current 
without major consequences to our bistability discussion. 


IV. MODEL PARAMETERS 


In order to solve Eq. we need to specify the matrix 
elements Tj in Eq.(|^. A simple and explicit model for 
qj is provided by Ref. [17]: 


Tj and the wavevectors 




1 


T • = t 


qj = ^(sin(^j,-cos(^j), 


( 12 ) 

with (fA = O 5 = 27r/3, ipc = 47 r/ 3 . This representa¬ 
tion is obtained for small twisting angles after perform¬ 
ing a 0 rotation of phase space in layer 2 (see details in 
Ref. [17]). It is also implicit in Eq.( 12 ) that the top and 
bottom graphene lattices have a common lattice point 
HZ]; a rigid horizontal translation between lattices adds 
an additional overall phase to the matrix Tj [18]. We 
stress, however, that relative phases in Tj do not alter 


in any significant way the physics of tunneling in Eq. 
Eurthermore, while the interlayer hopping amplitude t is 
sensitive to several parameters, e.g. twist angle m and 
the choice of dielectric material [20], its order of mag¬ 
nitude is mainly governed by the wavefunction overlap 
between the graphene layers. Such dependence will be 
described below within the WKB approximation. Equa¬ 
tions (§ and ( [T^ are expected to be accurate for twist 
angles 0 ^ 10°, and energies of 1 eV [19]. 

Eor an estimate below we use the value 0 = 2°. This 
defines an energy scale A = 0.37 eV. Eurthermore, we 
take a Lorentzian spectral function in Eq.(|6) for both 
layers, Ai^s(k,ci;) = 2r/ [{uj — shv-F\]^\)‘^ +1^ with the 
linewidth F ^ 10 meV. A finite linewidth F is necessary 
to have a finite value of the peak current when eVi 2 = A 
(see Fig|^. The temperature and Fermi level of the sys¬ 
tem were taken to be T = 0 and = 0, respectively. 
With reference to Eq. we define the interlayer con¬ 
ductance 


Gi2 = 5'^i 2, ^12 = 27rV 


|i| 


2 e2 


{hvY)‘^ h 


(13) 


where S is the surface area of the device. The value of gi 2 
is sensitive to the twist angle and the stacked dielectric 
material, if any, via the parameter t. Here we use gi 2 = 
Similar values of gi 2 were measured in 
resonant tunneling devices which contained four layers 
of BN in-between the graphene layers [12]. Eor Z, we 
consider a value of Z = G 12 /G 23 = 0.2. 


V. BISTABLE I-V CHARACTERISTICS 

The bistable I-V characteristics are shown in Eigj^. 
Eor a sufficiently large bias, eV ^ A, the current 
branches into two stable states. The low-resistance 
branch in Eigj^ corresponds to two layers at resonance 
(i.e., eVi 2 ~ A), whereas the high-resistance branch cor¬ 
responds to a non-resonant state (i.e., eVi 2 > A). We 
note that a third solution is also possible, indicated by a 
dashed line in the I-V response (see Eig{^). This solu¬ 
tion, however, is unstable given that a small perturbation 
in 6 V 12 will push the system away from this state. 

The bistable bias range can be estimated as SV ~ 
(^ 12 ^^ “^ 12 ^^) 7 ^ 23 , where /i 2 ^^ is the peak interlayer cur¬ 
rent and 1 ^ 2 ^^ is the valley interlayer current (see inset of 
Eigj^). To estimate /i 2 ^^ and Ii 2 ^\ we first note that 
the tunneling matrix element (k) varies, upon inte¬ 
gration in k-space, in the range 0 < \Tj^ (k)| < 2 t taking 
typical values |TJ^'(k)| t. Thus, it is a good approx¬ 
imation to take band and wavevector-independent tun¬ 
neling matrix elements (k)| = T. Eurthermore, in 
the typical case scenario the model parameters satisfy 
r(^ lOmeV) <C A(^ 0.1 — leV). With this in mind, 
the integration of Eq. allows an analytical expression 
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(b) 6^12 = A 




FIG. 3. Regions in k-space contributing to the resonant 
tunneling current for fixed V. These regions, indicated with 
black dashed lines at the intersection of the twisted Dirac 
cones, form conical surfaces in the k-plane: when eVi 2 < A 
the lines form hyperbolic curves, and when eVi 2 > A the 
lines form ellipsoidal curves. When eVi 2 = A, a van-Hove 
singularity in the tunneling density of states is obtained. As 
shown in (d), the non-resonant (high-resistance) bistable state 
(eVi 2 > A) can be Pauli-blocked by adjusting the doping. 
Doping thus affords a way to tune the current ratio between 
bistable branches in Fig[^. In this work it is assumed that 
the Dirac cones are aligned at V = 0, and that capacitance 
effects are neglected. Layers are labeled 1-3 as in Fig[^. 


to result in terms of line integrals in conical surfaces (see 
Figj^and the discussion in Appendix A). Using /ii ^2 = 0 
and Ui 2 = ^ 12 , we find that the non-resonant interlayer 
current takes the simple form 


luix) _ x^- 1/2 ,(vi) _ GuA 

/(y 4 e • 

Here x = eVi 2 /A > 1 and is the valley current 
obtained at a:; = •\/3/2. When eVyi/A = 1, however, the 
current is at resonance and reaches a maximum value 
which is sensitive to F. To leading order in F, we obtain 
(see Appendix B) 

(15) 


where F, in general, depends on the amount and type of 
disorder and/or temperature. Equations (14) and (15) 
yield e5Vj /S. 3v^T^Z[7rY^A/2F — l]/4. Importantly, 
very small values of Z (G 23 ^ G 12 ) make the bistable 
bias range negligibly small, whereas large values of Z 
(G 23 ^ G 12 ) would push the onset of the bistability re¬ 
gion to very large bias potentials. Optimally, values of 
Z ^ 1 and very small F would make the bistability effect 
more prominent. 

Achieving a large current ratio between bistable states 
is desirable for applications; this facilitates the reading 
process in a bistable device. From Eqs. in and ( p!^ , it is 


obtained that the current ratio between bistable branches 
is controlled by the parameter Z ^A/F. For realistic val¬ 
ues of disorder, this ratio can be in the 1-20 ballpark. It 
is interesting to note that these already high values can 
be boosted by means of Pauli blocking. As shown in 
Figsj^ and d, for sufficiently heavily doped samples, the 
non-resonant bistable state (but not the resonant one) is 
Pauli-blocked. The degree of the electrical current ratio 
enhancement depends on second order processes which 
assist tunneling, such as scattering with defects or dis¬ 
order. These second order processes are not considered 
here. 

The geometric control of Z, an appealing aspect of our 
system, can be understood from the Bardeen Transfer 
Hamiltonian Theory [22l [23]. In this theory, the inter¬ 
layer coupling t is calculated from the overlap of the 
wavefunctions of layers i and j in the barrier region, 
t = I2me) f dS • — '0jV'0*), with dS being 

a surface area element. Considering electrons tunnel¬ 
ing across a square potential barrier with a height much 
larger than the electron kinetic energy, a tunneling ma¬ 
trix element of the form t ex exp(—d^^/A) is obtained, 
where A is the WKB decay length defined above. The 
expression of Z in Eq. (§ results from assuming barriers 
between layers 1 and 2 and between layers 2 and 3 are of 
the same material, in combination with Eq.(13). 

Although electrostatic doping of the graphene layers 
is not essential for the physics that we describe, it is a 
convenient feature of bistability. In particular, for a fixed 
external bias potential, each bistable state exhibits differ¬ 
ent carrier concentrations. Thus, any in-plane measure¬ 
ment, such as conductance or magneto-transport, will 
be able to distinguish two distinct bistable states. In¬ 
deed, from the inset of Eig{^ we see that the interlayer 
bias potential for each bistable state differs by an amount 
SV 12 ^ A/e (see also discussion in the Appendix A). Tak¬ 
ing into account the capacitance of the layers, the induced 
carrier difference between both states is approximately 
Sn ^ K.A/A7Te^di2 (here the quantum capacitance is not 
included). Using 0 = 2°, k = and di 2 = 1 nm, we ob¬ 
tain a carrier density difference Sn ^ 10^^ cm“^ between 
stable states. These large carrier density differentials can 
be used as a smoking gun of intrinsic bistability. 


VI. OTHER GRAPHENE-BASED BISTABLE 
SYSTEMS 

Although we considered here for simplicity a two-step 
sequential tunneling structure where only one pair of lay¬ 
ers can be resonant, similar ideas apply to more com¬ 
plex structures. Interesting examples include a two-step 
resonant-resonant structure, opening the possibility for 
tristability or multi-step “cascade” devices. 

Einally, we also expect bistable I-V characteristics 
in twisted graphene trilayers in the absence of any di- 
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electric material. Indeed, incommensurability between 
graphene lattices already suppresses interlayer hybridiza¬ 
tion, regardless of the layers being spatially separated by 
a fraction of a nanometer, thus enabling the sequential 
tunneling regime [18]. Furthermore, the massless Dirac 
spectrum, and thus Eq. <© and the subsequent transport 
model, remain valid but with a modified Fermi velocity 
m- We stress, however, that stacked dielectric materials 
have two important advantages: (i) they enable tuning 
the interlayer coupling and (ii) they facilitate the inter¬ 
layer potential build-up in order to achieve a resonant 
behavior. 


VII. SUMMARY 

In summary, graphene-based van der Waals het¬ 
erostructures afford a new platform to realize devices 
with tunable I-V characteristics, in particular those with 
intrinsically bistable and hysteretic behavior. System 
parameters required to realize the bistable behavior are 
readily accessible in on-going experiments. The atomic 
scale interlayer distances can result in a fast response 
and large packing-densities, making these heterostruc¬ 
tures appealing for a variety of applications. 
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Appendix A: Capacitance effects 

In the main text, we fixed the Fermi energy /i^ of the 
different graphene layers to some constant value. A more 
refined model of the I-V response should, however, in¬ 
clude quantum capacitance effects so that Fermi energy 
is allowed to vary with V. Although the features of bista¬ 
bility are not significantly modified by such corrections, 
as shown below, carrier density differentials between the 
bistable states are a smoking gun of intrinsic bistability. 
These electrostatic considerations are discussed next. 

Here we numerically solve Eqs.([^-([TT|, assuming a 
thin device separated by dielectric barriers of thickness 
di 2 = (^23 = 1.4 nm (e.g. four layers of hBN) and di¬ 
electric constant k = 5. The procedure to solve the I-V 
response self-consistently is shown in Fig |Al^ , where rii 
and 77-2 are taken as independent variables [ns is obtained 
from Eq.([TQ|)], and 51 = I 12 —I23 in Er- (H is numerically 
calculated (color map). For fixed U, indicated with dot¬ 
ted isolines in Fig |Al^ , the self-consistent solutions to 
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FIG. Al. Self-consistent bistable solutions including quantum 
capacitance effects [Eqs.( 10)-(11)]. For fixed V, we find ni 
and 712 such that SI = I 12 — I 23 = 0. The bias isolines from 
Eq.([TT| are marked with dashed (V = 0) and dotted (finite 
V) lines, with an arrow pointing towards increasing V. The 
self-consistent I-V curve, obtained from the intersection of 
SI = 0 and the E-isolines in (a), is plotted in panel (b). 


the equilibrium equations are given by the pair ( 711 , 712 ) 
such that SI = 0. 

The resulting I-V response is shown in Fig |Al| 3. Im¬ 
portantly, the I-V characteristics are qualitatively simi¬ 
lar to those obtained in the main text by neglecting quan¬ 
tum capacitance effects. Furthermore, by inspection of 
the 77i and 77.2 axes in Fig |Al^ , we see that the differ¬ 
ence in carrier concentration Sn between each bistable 
state is on the order of Sn ^ cm“^. These car¬ 

rier concentration differences can easily be detected by 
lateral transport measurements and may act as clear fin¬ 
gerprints of intrinsic bistability. 


Appendix B: Analytic expressions for the peak and 
valley resonant tunneling current 


We derive here Eqs. 0 and ( [T^ of the main text, ob¬ 
tained under the assumption that the tunneling matrix 
elements Tj^' in Eq.(|^ are independent of the wavevector 
and band index, i.e., (k)| = T. Under this assump¬ 

tion, 1 12 depends only on the modulus of but not on its 
direction, and \Tj^ (k)P = 3T^. Given that F ^ A, 
when e^i 2 > A (non-resonant state) we can set F —0 
and thus take A^^ 5 (k, cj) = 27rS{uj — 5fiUF|k|). The two 
(5-functions appearing in Eq.(|^ can then be integrated in 
k-cj space, resulting in a one-dimensional integral along 


the contour of an ellipse: 

/ dk rh j 

(2^ J ^'^(‘^-s^^F|k|)(5(w-s'ai;F|k + q|) 




/ 




.A2 


sm 


167r3(?iVF)2 y(e$i2)2 - A2 

(Al) 

Here we denote a; by cD = cj + e^i 2 . In addition, the 
limits of integration on uj are given by uji = e^i 2 + /i 2 
and UJ 2 = /ii, whereas the limits of integration on (j) are 


7r/2, Xi> 1 

—l<Xi<l 

—7r/2, Xi < —1 


= 


2/ii^2 e^i2 


(A2) 

In obtaining Eq.( |Al[ ), we parametrized k-space using co¬ 
ordinates kx = kr sin (p/2 and ky = y^k‘^ — q‘^ cos (/)/2, 
with q conveniently aligned in the x-direction. The in¬ 
tegration over kr absorbs the first S function, setting 
kr = e^i 2 /hvF. Integration over uj absorbs the sec¬ 
ond S function, fixing the limits of integration ^ 1^2 in 
Eq.(A2). Importantly, because ^12 > A, the two S- 


functions in Eq.( Al) can only be non-zero simultaneously 


when s = — and s' = + (i.e. holes of layer 1 tunnel into 
electronic states of layer 2, see Fig|^. Using /ii ^2 = 0 
and V 12 = ^ 12 , Eqs.( |Al[ ) and ( |A2| ) result in Eq.([l^. 

When e^i 2 = A, it is necessary to restore the finite 
linewidth to the Lorentzian spectral function Ai^s(k, uj) = 
2F/ [(cj — shv-F\^\)‘^ + r^]. In this case, the integral for 
the tunneling current yields 


y f—T 

2 




^Ai,s(k,w)A2,s 

ZTT 


(k + q, w) = 


{hvF)^VTA 


/ duj \uj{uj — y 0{T/A) 

d U 2 


(A3) 

In obtaining Eq.( |A3| ), we transformed the integral of the 
spectral functions into a dimensionless integral of the 
form/res(e) = J (i^x{[/(x)^ + e] [ 5 f(x)^ + e]}“^ The func¬ 
tions / and g satisfy /(O) = 5^(0) = 0 and have a null 
Jacobian det[5x/, 9x^](0) = 0 (here e = F/A). It can 
be shown that Ues oc when e <C 1. An expansion 

to leading ord er in powers of e gives Eq.( |A3| ). Setting 

in Eq.(15) of 


/ii ^2 = 0 in Eq.(|A3|), the peak current 1 ^ 2 ^^ 


the main text is obtained. 
































